/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TANucleus
  \brief Info on a nucleus, but currently concentrated on 1N removal info, i.e.
  separation energies, energy levels, and c2s.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/05/25
  \date Last modified: 2022/05/25 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <iostream>
#include "TANucleus.h"
#include "TARs.h"
#include "TAMOMDIS.h"
#include "TAParaManager.h"
#include "TAException.h"
#include "TAYaml.h"

using std::cout;
using std::endl;

TANucleus::TANucleus(int a, double sp, double sn, double cs, double dcs) :
fA(a), fSp(sp), fSn(sn), fCSExp(cs), fDCSExp(dcs), fEk(-9999.), fPauli(-1),
fFNNParOpt(-1){}

TANucleus::TANucleus(const string &f){ // f: the yaml config file
  SetBaseDir(f.substr(0, f.find_last_of('/')));
  const TAYaml v(f);
  fA = v["nu"][0].as<double>();
  fSp = -9999.; // so as to distinguish between momdis and c.s., for subclass TANuMom
  if(v["nu"].size() == 5){
    fSp = v["nu"][1].as<double>();
    fSn = v["nu"][2].as<double>();
    fCSExp = v["nu"][3].as<double>();
    fDCSExp = v["nu"][4].as<double>();
  } // end if

  for(const auto &t : v["ex"]) fEx.push_back(t.as<double>());
  for(const auto &t : v["c2s"]) fC2S.push_back(t.as<double>());
  for(const auto &t : v["ob"]) fOrbit.push_back(t.as<string>());

  if(v["ek"]) fEk = v["ek"].as<double>();
  if(v["isPauli"]) fPauli = v["isPauli"].as<bool>();
  if(v["FNNParOpt"]) fPauli = v["FNNParOpt"].as<double>();
} // end ctor

// set the valence configurations
void TANucleus::SetOrbit(const initializer_list<const char *> &l){
  fOrbit.clear();
  for(const auto &t : l) fOrbit.push_back(t);
} // end member function SetOrbit

// set the excitation energies
void TANucleus::SetEx(const initializer_list<double> &l){
  fEx.clear();
  for(const auto &t : l) fEx.push_back(t);
} // end member function SetEx

// set the spectroscopic factors
void TANucleus::SetC2S(const initializer_list<double> &l){
  fC2S.clear();
  for(const auto &t : l) fC2S.push_back(t);
} // end member function SetC2S

void SetCS(const initializer_list<double> &l); ///< use external calculated cs
// use external calculated cs
void TANucleus::SetCS(const initializer_list<double> &l){
  fCS.clear();
  for(const auto &t : l) fCS.push_back(t);
} // end member function SetCS

const string &TANucleus::GetOrbit(int i) const{
  if(i < 0 || i >= int(fOrbit.size()))
    TAException::Error("TANucleus", "GetOrbit: index out of bounds");
  return fOrbit[i];
} // end member function GetOrbit

double TANucleus::GetEx(int i) const{
  if(i < 0 || i >= int(fEx.size()))
    TAException::Error("TANucleus", "GetEx: index out of bounds");
  return fEx[i];
} // end member function GetEx

double TANucleus::GetC2S(int i) const{
  if(i < 0 || i >= int(fC2S.size()))
    TAException::Error("TANucleus", "GetC2S: index out of bounds");
  return fC2S[i];
} // end member function GetC2S

void TANucleus::SetEk(double ek){
  if(ek <= 0.) TAException::Error("TANucleus", "SetEk: invalid Ek: %f", ek);
  fEk = ek;
} // end member function SetEk
void TANucleus::SetPauli(bool pauli){
  fPauli = pauli;
} // end member function SetPauli
void TANucleus::SetFNNParOpt(int opt){
  if(0 != opt || 1 != opt)
    TAException::Error("TANucleus", "SetFNNParOpt: invalid input: %d; 0: horiuchi; 1: LenziRay", opt);
  fFNNParOpt = opt;
} // end member function SetFNNParOpt

// calcualte Rs
void TANucleus::Rs(){
  if(-9999. == fSp) TAException::Error("TANucleus", "Rs: Sp not assigned!");
  const int n = fOrbit.size(), n1 = fEx.size(), n2 = fEx.size();
  if(n != n1 || n != n2) TAException::Error("TANucleus", "Rs: inconsistent vector sizes!");

  PARA->SetBaseDir(fBaseDir);
  // PARA->Clear(); // clean the previous text output of the base directory
  TARs rs(fA, fSp, fSn, fCSExp, fDCSExp); // the Rs calculator; true: -p; false: -n
  const bool exCS = int(fCS.size()) >= n; // whether use external CS or not
  for(int i = 0; i < n; i++){ // loop over all the valence configurations
    // set global parameters   fFNNParOpt: 0-> Horiuchi; 1-> LenziRay //
    if(0 == fFNNParOpt || 1 == fFNNParOpt) PARA->SetValue("FNNParOpt", fFNNParOpt);
    if(-9999. != fEk) PARA->SetValue("Ek", fEk); // MeV/u
    if(0 == fPauli || 1 == fPauli) PARA->SetValue("isPauli", bool(fPauli)); // 0: off; 1: on
    // do the calculation //
    char cfg[128]{}; rs.YamlName(fOrbit[i], cfg); // build the file name of the yaml
    // skip the calculation of Bound and S-Matrices if exCS is true
    TAMOMDIS mom(cfg, false, false, false, exCS); // do the calculation
    if(!i) cout << mom.banner() << endl; // print the reaction info
    rs.SetIsProtonRemoval(mom.IsProtonRemoval()); // info extracted from the yaml file
    // cout << "rs.IsProtonRemoval(): " << rs.IsProtonRemoval() << endl; getchar(); // DEBUG
    const double cs = exCS ? fCS[i] : mom.Total();
    rs.AddConfig(fOrbit[i], fEx[i], fC2S[i], cs); // register the configuration
  } // end for over i
  rs.Rs(); // calculate Rs
} // end member function Rs
